## This file is for analysing spatial paneldata set "numbered_tx_county_emp_wells" using R package "splm"

## Export the following output to a txt file.
sink(file = "~/HOLA/research/green-jobs/Rfile/creatJob_yearly.txt", append = FALSE, type = c("output", "message"), split = FALSE)

## clear the memory
rm(list = ls())

## Load library packages
library("spdep");         # Define listw class, and transform it into a sparse matrix
library("plm")            # Panel data model
library("splm");          # Spatial panel data model
library("xtable");        # Export results to Latex
library("foreign");       # Inport data from Stata,etc.
library("calibrate")      # Show labels to specific points in plots
library("ggplot2")


## Import data as a data.frame.
txData <- read.csv("~/HOLA/research/green-jobs/data/txData_yearly.csv", header = T)

## Look over the data set
names(txData);
dim(txData);
head(txData);
tail(txData);

## Creat identifier variable windcty/shalecty

txData$postwind <- as.numeric(txData$cumucap>50)
pwmat = matrix(txData$postwind, 254,11,byrow = T)
pwmat[rowSums(pwmat)>0] = 1
txData$windcty = c(t(pwmat))

txData$postshale <- as.numeric(txData$wells>15) 
psmat = matrix(txData$postshale, 254,11,byrow = T)
psmat[rowSums(psmat)>0] = 1
txData$shalecty = c(t(psmat))


## Create panel data frame
txPanel <- pdata.frame(txData, index = c("county","year"), drop.index = TRUE, row.names = TRUE)


## Examine the dist. of the data
summary(txPanel);
## stem(wells);
## dev.new()
## plot(sort(wells),pch=".")

## pdf("~/HOLA/research/green-jobs/latex/figures/hist_wells.pdf", height=6, width=8)
## hist(wells, prob = TRUE, col = "blue", main = NULL);
## lines(density(wells));
## rug(wells);
## dev.off()

## pdf("~/HOLA/research/green-jobs/latex/figures/dataPair.pdf", height=6, width=8)
## pairs(txPanel)
## dev.off()

## summary indicates the total variation of the variable and the share of this variation that is due to the individual and the time dimensions.
summary(txPanel$emp)
summary(txPanel$lemp)
summary(txPanel$wageR)
summary(txPanel$wageN)
summary(txPanel$wells)
summary(txPanel$cumucap)
## as.matrix indicates the total variation of the variable and the share of this variation that is due to the individual and the time dimensions. 
## head(lag(txPanel$emp, 0:2))
## head(diff(txPanel$emp), 15)
head(Within(txPanel$emp),15)
head(between(txPanel$emp),15)

##########################
## Impact on employment ##
##########################

## Basic regression model
## Dependent Variable: unemployment
## fm.unemp <- unempr ~ - log(wageR) - log(emp) - lag(log(wageR),1) + wells - lag(wells,1) + cumucap.pre - lag(cumucap.pre,1:2) - newCap - windInd - lag(windInd,1)
## fm.gmmunemp <- unempr ~ lag(unempr,1) - lag(wageR,1) + wells - lag(wells,1) + cumucap.pre - lag(cumucap.pre,1)- windInd - lag(windInd,1)|lag(unempr,2:11)
txData$date = as.yearmon(txData$year)
xyplot(wells~date|county,data = txData[txData$county == "Harris",])

fm.emp <- emp ~ cumuwells
fm.logemp <- log(emp) ~ wells + cumucap.pre
fm.gmmemp <- log(emp) ~ lag(log(emp),1) + wells + cumucap.pre|lag(log(emp),2:11)

## Pooling effects:
job.pool <- plm(fm.logemp,txPanel,model = "pooling")
## Fixed effects:
## One-way
job.fe <- plm(fm.logemp,txPanel,model = "within")
job.re <- plm(fm.logemp,txPanel,model = "random")
## Two-way
job.twfe <- plm(fm.emp,txPanel,model = "within", effect = "twoways")
summary(job.twfe);
job.twre <- plm(fm.emp,txPanel,model = "random",random.method = "amemiya", effect = "twoways")
summary(job.twre)

fe.spatial = fixef(job.twfe, type = "dmean");

bmp("~/HOLA/research/green-jobs/latex/figures/spatialEf.bmp", height=6, width=8, units = "in",res = 300)
plot(fe.spatial,main = "County-specific fixed effects", xlab = "county #", ylab = "Fixed effects, # of employment")
textxy(fe.spatial, names(fe.spatial))
dev.off()

pdf("~/HOLA/research/green-jobs/latex/figures/spatialEf.pdf", height=6, width=8)
plot(fe.spatial,main = "County-specific fixed effects", xlab = "county #", ylab = "Fixed effects, # of employment")
textxy(fe.spatial, names(fe.spatial))
dev.off()

fe.time = fixef(job.twfe, effect = "time");

bmp("~/HOLA/research/green-jobs/latex/figures/timeEf.bmp", height=6, width=8, units = "in",res = 300)
plot(fe.time,main = "Time effects", xlab = "Year", ylab = "Time effects, , # of employment", type = 'o')
dev.off()


pdf("~/HOLA/research/green-jobs/latex/figures/timeEf.pdf", height=6, width=8)
plot(fe.time,main = "Time effects", xlab = "Year", ylab = "Time effects, , # of employment", type = 'o')
dev.off()

job.pv <- pvcm(fm.logemp,txPanel,model = "within")
summary(job.pv)
pooltest(job.fe,job.pv,effect = T)
plmtest(job.fe,effect = "time")
pFtest(job.twfe,job.pool)

jobloss.pool <- plm(fm.unemp,txPanel,model = "pooling")
## Fixed effects:
## One-way
jobless.fe <- plm(fm.unemp,txPanel,model = "within")
jobless.re <- plm(fm.unemp,txPanel,model = "random")
## Two-way
jobless.twfe <- plm(fm.unemp,txPanel,model = "within", effect = "twoways")
summary(jobless.twfe);
jobless.twre <- plm(fm.unemp,txPanel,model = "random",random.method = "amemiya", effect = "twoways")
summary(jobless.twre)



job.logtwts <- pgmm(formula = fm.gmmemp,data = txPanel, model = "twosteps", effect = "twoways")
summary(job.logtwts);
jobloss.logtwts <- pgmm(formula = fm.gmmunemp,data = txPanel, model = "twosteps", effect = "twoways")
summary(jobloss.logtwts);

## DID Method

## Import spatial weights matrix.
wmat = as.matrix(read.csv("~/HOLA/research/green-jobs/data/wmat/adjmatrix.csv", header = TRUE, row.names = 1))


wellmat = as.matrix(txPanel$wells)
windmat = as.matrix(txPanel$cumucap)
pairshale <- matrix(nrow = 100,ncol = 4)
pairwind <- matrix(nrow = 200,ncol = 4)
k = 1;
l = 1;
for (i in 1:254) {
  for (j in 1:254) {
    if (wmat[i,j] > 0 && any(wellmat[i,] > 10) && all(wellmat[j,] == 0)) {
        pairshale[k,] = c(index(wmat)[i],rownames(wmat)[i], index(wmat)[j], colnames(wmat)[j])
        k = k+1 
      }
    if (wmat[i,j] > 0 && any(windmat[i,] > 50) && all(windmat[j,] == 0)) {
        pairwind[l,] = c(index(wmat)[i],rownames(wmat)[i], index(wmat)[j], colnames(wmat)[j])
        l = l+1 
      }
  }
}

pairshale = na.omit(pairshale)
pairwind = na.omit(pairwind)

pair <- pairshale;


fm.shaledid <- lemp ~  shalecty + postshale + shalecty*postshale
fm.winddid <- lemp ~  windcty + postwind + windcty*postwind
coef.did <- matrix(nrow = dim(pair)[1],ncol = 2)

k = 1;
for (i in 1:dim(pair)[1]) {
txPair = subset(txPanel, countyID == pair[i] | countyID == pair[i,3])
## txPair$postwind[txPair$countyID == pair[i,3]] = txPair$postwind[txPair$countyID == pair[i]]
## emp.did  = lm(fm.winddid, txPair)
txPair$postshale[txPair$countyID == pair[i,3]] = txPair$postshale[txPair$countyID == pair[i]]
emp.did  = lm(fm.shaledid, txPair)

summary(emp.did)
coef.did[k,1] = coef(emp.did)[4]
if (!is.na(coef.did[k,1])) {
## coef.did[k,2] = summary(emp.did)$coefficients["windcty:postwind","Pr(>|t|)"]
coef.did[k,2] = summary(emp.did)$coefficients["shalecty:postshale","Pr(>|t|)"]
}
k = k+1
}
coef.did = na.omit(coef.did)
coef.did.pvalue = coef.did[coef.did[,2]<0.001,]
coefDID = mean(coef.did[,1])
coefDIDp = mean(coef.did.pvalue[,1])
txPair = subset(txData, countyID == pair[i] | countyID == pair[i,3])
xyplot(unempr~year, data = txPair, groups = countyID, type = "o", auto.key=list(title="County", space = "right", cex=1.0))


txPair$wind <- as.numeric(txPair$cumucap.pre > 0);
txData$lemp <- log(txData$emp)

a = sapply(subset(txPair, shalecty == 0 & postshale == 0, select = lemp), mean)
b = sapply(subset(txPair, shalecty > 0 & postshale == 0, select = lemp), mean)
c = sapply(subset(txPair, shalecty == 0 & postshale > 0, select = lemp), mean)
d = sapply(subset(txPair, shalecty >0  & postshale > 0, select = lemp), mean)

d-c-(b-a)

e = sapply(subset(txPair, post06 ==0 & wind == 0, select = lemp), mean)
f = sapply(subset(txPair, post06 ==0 & wind == 1, select = log(emp)), mean)
g = sapply(subset(txPair, post06 ==1 & wind == 0, select = log(emp)), mean)
h = sapply(subset(txPair, post06 ==1 & wind == 1, select = log(emp)), mean)

h-g-(f-e)

####################
## Impact on wage ##
####################
## Basic regression models
fm.wage <- wageR ~ -lag(wageR,1) - emp + wells - lag(wells,1) + cumucap.pre - lag(cumucap.pre,1) - windInd - lag(windInd,1) - newCap - lag(newCap,1) 
fm.logwage <- log(wageR) ~ - lag(log(wageR),1) - log(emp) + wells - lag(wells,1) + cumucap.pre - lag(cumucap.pre,1) - windInd - lag(windInd,1)
fm.gmmwage <- log(wageR) ~ + lag(log(wageR),1) - lag(log(emp)) + wells + lag(wells,1) + cumucap.pre + lag(cumucap.pre,1) - windInd - lag(windInd,1)|lag(log(wageR),2:11)

## Pooling effects:
pay.pool <- plm(fm.wage,txPanel,model = "pooling")
## Fixed effects:
## One-way
pay.fe <- plm(fm.wage,txPanel,model = "within")
pay.re <- plm(fm.wage,txPanel,model = "random")
## Two-ways
pay.twfe <- plm(fm.wage,txPanel,model = "within", effect = "twoways")
summary(pay.twfe);
pay.twre <- plm(fm.wage,txPanel,model = "random",random.method = "amemiya", effect = "twoways")
summary(pay.twre)

pay.pv <- pvcm(fm.wage,txPanel,model = "within")
summary(pay.pv)
pooltest(pay.fe,pay.pv,effect = T)
plmtest(pay.fe,effect = "time")
pFtest(pay.twfe,pay.pool)

pay.logtwfe <- plm(fm.logwage,txPanel,model = "within", effect = "twoways")
summary(pay.logtwfe);

pay.logtwts <- pgmm(formula = fm.gmmwage,data = txPanel, model = "twosteps", effect = "twoways")
summary(pay.logtwts);
## Random effects:

## One-way effect, from Swamy and Arora (1972)
job.re <- plm(fm.logemp,txPanel,model = "random")
summary(job.re)
## Two-way effects, from Amemiya (1971),
job.twre <- plm(fm.logemp,txPanel,model = "random",random.method = "amemiya", effect = "twoways")
summary(job.twre)
attach(job.twre)
pdf("~/HOLA/research/green-jobs/figures/twre.pdf", height=6, width=6)
plot(txPanel$wells,txPanel$emp, main = "Random effects model")
abline(coefficients, lwd = 2, col = "blue")
mtext(bquote( emp == .(coefficients[2]) * txPanel$wells + .(coefficients[1])), side=3, line=0) 
dev.off()
detach()

## Remove two outliers : Harris and Dallas
txDataR <- subset(txData, county != "Harris" & county != "Dallas")
txPanelR <- pdata.frame(txDataR, index = c("county", "year"), drop.index = TRUE, row.names = TRUE)

job.retwfe <- plm(fm.logemp,txPanelR, model = "within", effect = "twoways")
summary(job.retwfe);

job.retwre <- plm(fm.logemp,txPanelR,model = "random",random.method = "amemiya", effect = "twoways")
summary(job.retwre)

## For counties which have well drilling activities, analyse the impact of gas, oil and other wells.

txShaleWell <- subset(txData, wells != 0)
txShaleWell <- pdata.frame(txShaleWell, index = c("county", "year"), drop.index = TRUE, row.names = TRUE)
fm.emp.welltype <- emp ~ gas + oil

well.twfe <- plm(fm.emp.welltype,txShaleWell, model = "within", effect = "twoways")
summary(well.twfe);

well.re <- plm(fm.emp.welltype,txShaleWell, model = "random")
summary(well.re);


txShaleCty = txData
for (i in 1:254) {
  if (all(wellmat[i,] == 0)) {
   cty =  dimnames(wellmat) [[1]][i]
  txShaleCty = subset(txShaleCty, txShaleCty$county != cty)
 }
}
txShaleCty <- pdata.frame(txShaleCty, index = c("county", "year"), drop.index = TRUE, row.names = TRUE)

cty.twfe <- plm(fm.emp.welltype,txShaleCty, model = "within", effect = "twoways")
summary(cty.twfe);

cty.twre <- plm(fm.emp.welltype,txShaleCty, model = "random", effect = "twoways")
summary(cty.twre);

wt.twfe <- plm(fm.emp.welltype,txPanel, model = "within", effect = "twoways")
summary(wt.twfe);

wt.twre <- plm(fm.emp.welltype,txPanel, model = "random", effect = "twoways")
summary(wt.twre);

phtest(cty.twfe,cty.twre)
phtest(wt.twfe,wt.twre)
## Tests
## Tests for individual and time effects
plmtest(job.pool,effect = "twoways", type = "bp")
plmtest(job.pool,effect = "individual",type = "bp")
plmtest(job.pool,effect = "time",type = "bp")
pFtest(job.fe,job.pool)
pFtest(job.twfe,job.pool)

## Hausman test
phtest(job.fe,job.re)
phtest(job.twfe,job.twre)
phtest(job.retwfe,job.retwre)

## Tests of serial correlation
## joint/conditional LM test, test for ramdom effects and serial correlation under normality and homoskedasticity of the idiosyncratic errors
pbsytest(fm.logemp, txPanel, test = "j")
pbsytest(fm.logemp, txPanel, test = "re")
pbsytest(fm.logemp, txPanel)
## General serial correlation tests
pbltest(fm.logemp, txData, alternative = "onesided")
pbgtest(job.twfe, order = 2)
## For "short" panels with small T and large n
pwartest(fm.logemp, txData)
pwfdtest(fm.logemp, txData)
## Test for cross-sectional dependence
pcdtest(fm.logemp, txData, model = "within", effect = "twoways")
## Robust cov matrix estimation
coeftest(job.twre,vcovHC)


#######################
## ML Implementation ##
#######################



## spacial weights matrix
wmatl <- mat2listw(wmat, style = "W");

## Spatial error model, random effects
## Results for the Kapoor et al. (2007)
errorre <- spml(formula = fm.wage, data = txData, index = c("county","year"),listw = wmatl,  model = "random", lag = F, spatial.error = "kkp")
summary(errorre)

lagre <- spml(formula = fm.wage, data = txData, index = c("county","year"), listw = wmatl,  model = "random", lag = T, spatial.error = "none")
summary(lagre)

## Spatial lag & error model, random effects
errlagre <- spml(formula = fm.wage, data = txData, index = c("county","year"), listw = wmatl,  model = "random", lag = TRUE, spatial.error = "kkp")
summary(errlagre)

## Spatial individual effects mu_i only, fixed effects
spacefe <- spml(formula = fm.logemp, data = txData, index = c("county","year"), listw = wmatl, model = "within", effect = "individual", method = "eigen", na.action = na.fail, quiet = TRUE, zero.policy = NULL, interval = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE)
summary(spacefe)

## Time effects gamma_t only, fixed effects
timefe <- spml(formula = fm.logemp, data = txData, index = c("county","year"), listw = wmatl, model = "within", effect = "time", method = "eigen", na.action = na.fail, quiet = TRUE, zero.policy = NULL, interval = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE)
summary(timefe)

## Two-way effects, fixed effects
lagtwfe <- spml(formula = fm.logemp, data = txData, index = c("county","year"), listw = wmatl, model = "within", lag = T, spatial.error = "none", effect = "twoways", method = "eigen", na.action = na.fail, quiet = TRUE, zero.policy = NULL, interval = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE)
summary(lagtwfe)

errtwfe <- spml(formula = fm.logemp, data = txData, index = c("county","year"), listw = wmatl, model = "within", lag = F, spatial.error = "b", effect = "twoways", method = "eigen", na.action = na.fail, quiet = TRUE, zero.policy = NULL, interval = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE)
summary(errtwfe)

## Two-way spatial lag & error model, fixed effects
errlagtwfe <- spml(formula = fm.logemp, data = txData, index = c("county","year"), listw = wmatl, model = "within", lag = T, spatial.error = "b", effect = "twoways", method = "eigen", na.action = na.fail, quiet = TRUE, zero.policy = NULL, interval = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE)

summary(errlagtwfe)

eff <- effects(errtwfe)
eff

## #######################
## ## GM Implementation ##
## #######################

## ## Spatial error model, random effects

GM.error <- spgm(formula = fm.logemp, data = txData, index = c("county","year"), listw = wmat, moments = "initial", model = "random", spatial.error = TRUE)

summary(GM.error)

##Spatial error & lag  model, random effects
GM.full <- spgm(formula = fm.logemp, data = txData, index = c("county","year"), listw = wmat, lag = TRUE, moments = "fullweights", model = "random", spatial.error = TRUE)

summary(GM.full)

## #Spatial error model, fixed effects#
## GM_error <- spgm(formula = fm.emp, data = txData, lag = TRUE, listw = wmat, model = "within", spatial.error = TRUE)

## summary(GM_error)

#Testing#

## LM test, no random effects assuming no spatial correlation
## H0: sigma^2_mu = 0|lambda = 0
test1 <- bsktest(x = fm.logemp, data = txData, index = c("county","year"),listw = wmatl, test = "LM1")
print(class(test1))
test1

## LM test, no spatial correlation assuming no random effects
## H0: lambda = 0|sigma^2_mu = 0
test2 <- bsktest(x = fm.wage, data = txData, index = c("county","year"), listw = wmatl, test = "LM2")
test2

## Joint LM test
test3 <- bsktest(x = fm.logemp, data = txData,index = c("county","year"), listw = wmatl, test = "LMH")
test3

## Conditional LM_lambda test#
## H0: lambda = 0| sigma^2_mu >= 0 
bsktest(x = fm.logemp, data = txData, index = c("county","year"), listw = wmatl, test = "CLMlambda");

## Conditional LM_mu test#
## H0: sigma^2_mu = 0| lambda may or may not be zero
bsktest(x = fm.wage, data = txData, index = c("county","year"), listw = wmatl, test = "CLMmu");


## #Spatial Hausman test#
## testH <- sphtest(x = fm.wage, data = txData, index = c("county","year"),listw = wmatl, spatial.model = "error", method = "GM")
## ## testH

mod1 <- spgm(formula = fm.wage, data = txData, index = c("county","year"), listw = wmat, lag = T, moments = "fullweights", model = "random", spatial.error = T)

mod2 <- spgm(formula = fm.wage, data = txData, index = c("county","year"), listw = wmat, lag = T, moments = "fullweights", model = "within", spatial.error = T)

test2 <- sphtest(x = mod1, x2 = mod2)
test2

## #Linear hypothesis testing# 
## coeftest(errlagtwofe)

sink()









data(oldcol)
     Vx <- lag.listw(nb2listw(COL.nb, style="W"), COL.OLD$CRIME)
     plot(Vx, COL.OLD$CRIME)
     plot(ecdf(COL.OLD$CRIME))
     plot(ecdf(Vx), add=TRUE, col.points="red", col.hor="red")
    is.na(COL.OLD$CRIME[5]) <- TRUE
     VxNA <- lag.listw(nb2listw(COL.nb, style="W"), COL.OLD$CRIME, NAOK=TRUE)
